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We study one-dimensional, interacting, gapped fermionic systems described by variants of the 
Peierls-Hubbard model, and characterize their phases via a topological invariant constructed out of 
their Green's functions. We demonstrate that the existence of topologically protected, zero-energy 
states at the boundaries of these systems can be tied to the value of the topological invariant, just like 
when working with the conventional, noninteracting topological insulators. We use a combination of 
analytical methods and the numerical density matrix renormalization group method to calculate the 
values of the topological invariant throughout the phase diagrams of these systems, thus deducing 
when topologically protected boundary states are present. We are also able to study topological 
states in spin systems because, deep in the Mott insulating regime, these fermionic systems reduce 
to spin chains. In this way, we associate the zero-energy states at the end of an antiferromagnetic 
spin-one Heisenberg chain with a topological invariant equal to 2. 
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I. INTRODUCTION 

Topological insulators — free fermionic systems with 
topological band structure 1-6 — are now very well under- 
stood. A band structure can be called topological if 
it has a nonvanishing topological invariant, the Chern 
number first proposed for the integer quantum Hall ef- 
fect being the simplest example. 7 These invariants imply 
zero-energy boundary states; it is these boundary states 
that distinguish the topological insulators from their non- 
topological counterparts and which are crucial for their 
physical properties. However, a number of topological 
states of matter have been discovered whose existence 
requires interactions. Prominent among these are the 
"topologically ordered" states in two dimensions as de- 
fined by X.-G. Wen, inspired by the fractional quantum 
Hall effect. 8 Another example is the Haldanc state of 
the spin-one antiferromagnetic Heisenberg chain. 9 Like 
noninteracting topological insulators, these states are 
bulk-incompressible with zero-energy excitations at the 
boundary. However, they are substantially different from 
the topological insulators in the need for interactions, 
and their excitations are often fractionalized relative to 
the underlying, microscopic degrees of freedom as, e.g., 
Laughlin's fractionally charged excitations in the frac- 
tional quantum Hall effect and the spin- 1/2 boundary 
states of the Haldane state of the spin-1 Heisenberg chain. 
It is therefore natural to ask whether there is a connec- 
tion between topological band structures and interacting 
topological states. 

Here, we explore one connection proposed in the lit- 
erature in recent years, 10-14 where one computes topo- 
logical invariants of the single-particle Green's function 
rather than the single-particle Hamiltonian. These in- 
variants coincide in the absence of interactions but, un- 
like single-particle Hamiltonians, single-particle Green's 
functions continue to exist even when interactions are 



present. Their topological invariants thus generalize the 
free invariants to generic, interacting systems. 

In this paper, we use this approach to compute explic- 
itly a topological invariant for one-dimensional, interact- 
ing, fermionic systems, and connect its value to the pres- 
ence or absence of topologically protected, zero-energy 
boundary states. In particular, we study Hubbard mod- 
els with dimerization, as well as spin chains, which can 
be understood as Mott-insulating phases of the Hubbard 
models in the strongly interacting regime. Where possi- 
ble, we calculate invariants analytically but, in the gen- 
eral case, we rely on the numerical density matrix renor- 
malization group method 16-18 and its time-dependent ex- 
tension (adaptive t-DMRG 19 ' 20 ) to compute the invariant 
as a function of the system parameters. Being an integer, 
the invariant is not very sensitive to the errors inherent 
in our numerical approach and so can be determined pre- 
cisely with moderate computational effort. The ultimate 
goal of this paper is not so much to calculate the phase 
diagram of the models we study as to illustrate the utility 
of the Green's-function method of topological invariants 
when applied to interacting, one-dimensional, fermionic 
systems. 

While we find this approach to be useful, some cau- 
tion is needed in the interpretation of these interacting 
topological invariants. First, in the absence of interac- 
tions the topological invariants frequently measure the 
linear response to external (frequently electromagnetic) 
perturbations. 7 ' 21 For example, the Chern number char- 
acterizing the integer quantum Hall states is proportional 
to the Hall conductance a xy of these systems. However, 
this connection is not guaranteed for the Green's func- 
tion invariant of an interacting system. Second, inter- 
actions introduce a novel possibility for the physics at 
the boundary not present in free systems. In particu- 
lar, in the free system, a nontrivial value for the bulk 
invariant implies zero-energy (single particle) excitations 
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at the boundary, which formally appear as poles of the 
single-particle Green's function. In an interacting sys- 
tem, a nontrivial bulk invariant is also consistent with 
zero-energy zeros of the Green's function, 14 in addition 
to poles, 10 ' 22 which indicates a complete loss of coherence 
of the single-particle degrees of freedom. We will see this 
explicitly in examples. 

In the bulk, this behavior — a singularity of the 
Green's function without a corresponding zero-energy 
single-particle state — means that the topological invari- 
ant can change its value (discontinuously) as parameters 
are varied without passing through a phase transition. 
This means that, unlike for free systems, these topolog- 
ical invariants do not necessarily correspond in a sim- 
ple way to phases of matter. Fortunately, recent works 
have addressed this issue for one-dimensional fermion 
systems. 23-25 These authors find that there are only four 
topologically distinct phases, characterized by the value 
of the topological invariant modulo 4. [Note that the 
authors of these papers also consider fermionic systems 
for which the number of particles is not conserved, i.e., 
without U(l) symmetry, which leads to 8 distinct phases, 
with Z§ structure. In this paper we restrict ourselves to 
systems with U(l) symmetry and Z4 structure.] Given 
this knowledge, our computations of topological invari- 
ants provide unambiguous determination of the topolog- 
ical phase for each model we consider. 

We apply these ideas first to the one-dimensional 
Peierls-Hubbard model of spin- 1/2 fermions with dimer- 
ized hopping and on-site Hubbard repulsion (see, e.g., 
Ref. 26). Its invariant can take on the values or 2. 
We also briefly treat a system of two coupled Hubbard 
chains, which realizes the Fidkowski-Kitaev model 23 and 
whose invariant can take on values which are multiples 
of four; i.e., this model has only one phase. We calculate 
the interacting invariants for these models using both an- 
alytical arguments and the t-DMRG method. We point 
out that the existence of the boundary states in these 
models can indeed be captured by these invariants, in ac- 
cordance with the bulk-boundary correspondence. 22 We 
further point out that when the invariant is a multiple 
of four, the disappearance of the boundary states is con- 
sistent with the bulk-boundary correspondence, thanks 
to the replacement of zero-energy boundary states by ze- 
ros in the Green's function. Finally, we analyze numer- 
ically and with some analytical arguments a variant of 
the fermionic Peierls-Hubbard chain whose parameters 
are adjusted so that the system is in a spin-one Haldane 
phase. We find that the interacting topological invari- 
ant is equal to 2, showing that the spin-one Heisenberg 
chain possesses the same boundary states as the Peierls- 
Hubbard model. 

The rest of the paper is organized as follows. In Sec. II, 
we introduce the topological invariant for interacting one- 
dimensional systems. In Sec. Ill, we introduce the one- 
dimensional Peierls-Hubbard model and discuss its phase 
diagram, topological invariant, and boundary states ana- 
lytically. In Sec. IV, we complete the study of the phase 



diagram of the Peierls-Hubbard model and the topologi- 
cal invariant using the DMRG method. In addition, we 
analyze the boundary states in a variant of the Peierls- 
Hubbard model with spin interactions that is equivalent 
to the spin-one Haldane chain. In Sec. V, we show that 
the absence of phase transitions in the Fidkowski-Kitaev 
model is compatible with a changing invariant precisely 
because the Green's function acquires zeros at a certain 
point of the phase diagram, leading to a change of the in- 
variant by a multiple of 4. Finally, in Sec. VI, we present 
our conclusions and outlook. 



II. ONE-DIMENSIONAL TOPOLOGICAL 
INSULATORS AND THEIR TOPOLOGICAL 
INVARIANT 

Consider a one-dimensional fermionic system. Two 
types of topological invariants are known to exist for 
these systems (in the absence of interactions), one with 
Z structure and one with Z2 structure. In this paper we 
concentrate on those systems whose invariant is an inte- 
ger Z. To allow for the existence of the topological invari- 
ant, the imaginary-time single-particle Green's function 
G(k,ui) of such systems must possess the symmetry 14 ' 15 

EG(fe,o;)S = -G(fc,-w), (1) 

where £ is some unitary matrix whose square is 1. Here uj 
is the imaginary frequency and k is the wave vector. This 
symmetry is usually referred to as a chiral symmetry. 1 ' 14 
In the absence of interactions, it occurs when a par- 
ticle moves on a bipartite lattice; in one dimension, 
this symmetry is present for a tight-binding model with 
nearest-neighbor hopping only. In the presence of inter- 
actions, it appears as a combination of particle-hole and 
time-reversal transformations and is also very ubiquitous. 
In particular, adding Hubbard-type interactions to the 
tight-binding model with sublattice symmetry preserves 
Eq. (1). 

A Green's function with the property Eq. (1) is char- 
acterized by a topological invariant. Defining 

ff(*)=G(fc,w)L =0 , (2) 

we can write the invariant as 12 ' 13,27 

^-^9- (3) 

Here the trace is taken over the matrix indices of g (which 
label the bands and spin indices of the model we study) . 
Written in this form, this invariant exists whether or not 
interactions are present. One comment is in order: this 
definition of the topological invariant differs from the one 
adopted in Ref. 23 by a factor of 2, since, as was pointed 
out in Sec. I, we work with complex (or Dirac) fermions 
with a conserved number of particles, while Ref. 23 works 
with real (or Majorana) fermions. The subscript 1 in Ni 
refers to the one-dimensional space in which it is defined. 
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It is straightforward to see that iVi is topological, that 
is, it does not change if one changes g(k) slightly. Indeed, 
if g depends on some parameter a (a coupling constant 
in the Hamiltonian, for example), the derivative dNi/da 
can be found to be 

i.e., as an integral over a total derivative. To show that 
Eq. (4) holds, one must take advantage of Eq. (1), which 
implies, together with Eq. (2), that g anticommutes with 
S, 

£.9 = -£<?• (5) 

It is well-known that a basis always exists for which the 
matrix E takes the form 

s = (J -i) • w 

where 1 stands for an identity matrix. This allows us 
to rewrite the topological invariant, Eq. (3), in a slightly 
different form by noting that thanks to the condition (5) 
as well as to Eq. (6), the matrix g must have the off- 
diagonal structure 

a-( ° v{k A (7) 

where v(k) is some generic matrix. Substituting Eq. (7) 
into Eq. (3), we find 

Ni=tr J ^d k logv(k)=Yl J ^dklog z n (k), 

(8) 

where z n (k) are the eigenvalues of v(k). We see that Ni 
simply counts the number of eigenvalues of v that wind 
around the origin of the complex plane as a function of 
k. 

In practical applications below, v(k) will often be just 
a number or a diagonal matrix, and z n (k) will be very 
straightforward to identify. 

In the absence of interactions, it is possible to relate 
G and g directly to the Hamiltonian. Indeed, a generic 
noninteracting Hamiltonian looks like 

H = Y,^cic^ (9) 

a/3 

where c and & are annihilation and creation operators, 
and the indices a, (3 refer to lattice sites, spin, and flavor 
of fermions, if any. Its Green's function is given by 

G= [iu-nr 1 , g = -U~ 1 . (10) 

With this identification, the invariant, Eq. (3), becomes 
the well-known one-dimensional version of the topological 
invariant for noninteracting fermionic systems with chiral 



symmetry. 1 It is used, for example, to identify topological 
phases and boundary states of fermionic chains such as 
those studied in Ref. 28 (whose boundary states are often 
referred to as Su-Schricffcr-Hceger solitons). 

Once interactions are turned on, however, simple ex- 
pressions such as Eq. (10) are no longer available. The 
utility of the topological invariant Eq. (3) lies in the fol- 
lowing. First, in the absence of interactions, the only 
way for iVi to change is if g becomes singular at some 
momentum k. This can happen only if the system has 
zero-energy excitations as follows from Eq. (10), thus im- 
plying a quantum phase transition. In the presence of 
interactions, N\ can also change if g acquires zero eigen- 
values at some k. (This is impossible in the absence of 
interactions as follows from Eq. (10) assuming that the 
Hamiltonian is bounded; see also Ref. 29, which relaxes 
this assumption.) This is the origin of the possibility that 
Ni may change value even in the absence of a quantum 
phase transition. 

Second, suppose we have two adjacent domains where 
the invariant Ni takes on two different values, in the 
right domain and in the left one. Then one can show 
that 22 

N?-N[ =tr£, (11) 

where the trace is evaluated in the Hilbert space of the 
chain with an open boundary. This constitutes the bulk- 
boundary correspondence for this type of topological in- 
sulator, and tr £ can be viewed as a boundary topological 
invariant, as we will see below. The derivation of Eq. (11) 
(see, in particular, Appendix B of Ref. 22) is based on 
the algebraic manipulations of the function g, indepen- 
dent of its physical meaning (in particular, independent 
of whether interactions are present). 

To see the implication of this relation for the boundary 
states, we calculate this trace in the basis of the eigen- 
states of g a p, the zero- frequency Green's function of the 
open chain (which is thus not translationally invariant 
and so g cannot be reduced to just a function of the 
momentum k). Every eigenstate ip n of g with nonzero 
eigenvalue A n , 

gipn=Klpn, (12) 

has a conjugate eigenstate with an opposite eigenvalue, 
gT,ip n = -X n T,tp n , (13) 
as a consequence of Eq. (5). It follows that 

C^n = (14) 

because ip n and T,ip n are both eigenstates of g with op- 
posite eigenvalues. Thus eigenstates with nonzero eigen- 
values A rl do not contribute to tr S. Only the eigenstates 
of g with A„ either infinite (poles of the Green's func- 
tion at to = 0) or zero (zeros of the Green's function at 
u) = 0) contribute to the trace. So we see that tr E must 
be counting zero-energy states present on the boundary 
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between two topological insulators, or possibly zeros of 
Green's functions if there are interactions, justifying its 
earlier characterization as a boundary topological invari- 
ant. 

Moreover, all such eigenstates with A„ either infinite 
or zero are also eigenstates of £ with eigenvalues that 
are either +1 or — 1. To see this, we observe that if ip n 
is an eigenstate with a zero eigenvalue, then so is £f/Vi- 
By forming linear combinations tp n ± T,ijj n , we construct 
the eigenstates that are also eigenstates of £ with the 
promised eigenvalues. Similar arguments can be given 
for eigenstates with infinite A„. Thus, tr£ counts these 
types of eigenstates with appropriate signs. We see that 
Eq. (11) simply tells us that a boundary between two 
topological insulators with different values of the topolog- 
ical invariant must have either some zero energy states, 
or some zeros of the Green's function, or perhaps both. 

We will see that in practical applications of Eq. (11) 
a third possibility can arise. Namely, the existence of 
zero-energy excitations signifies that the ground state is 
degenerate and the system must pick one of the degen- 
erate states spontaneously. The ground state chosen in 
this way may break the chiral symmetry and Eq. (1) may 
then break down. In this case, it may happen that the 
system has no single-particle zero-energy excitations and 
no zeros of the Green's function, since Eq. (11) is simply 
no longer valid. This happens in the dimerized Hubbard 
model studied in the next section under certain condi- 
tions. However, the system still has zero-energy bound- 
ary states, which is reflected in the existence of more 
than one ground state. It is just that the Green's func- 
tion, which is sensitive to the single-particle excitations 
only, does not see those zero-energy excitations corre- 
sponding to the multiple ground states, which are actu- 
ally of particle-hole type in this case. The main conclu- 
sion remains the same: under these conditions a nonzero 
bulk topological invariant implies zero-energy states at 
the edge. 

Finally, as we have pointed out in Sec. I, at the bound- 
ary of a topological system with the invariant iVi equal to 
modulo 4, there should be no zero-energy states. That 
means that, instead, there will only be zeros, in accor- 
dance with Eq. (11). The same is true for the boundary 
between two insulators with the difference of Ni across 
the boundary equal to modulo 4. The converse is also 
true: if the invariant changes by a number other than 
(mod 4) across a boundary, there will always be some 
zero-energy states. This follows from the arguments of 
Refs. 23-25 as well as for the reasons presented above. 
One interesting possibility is that the boundary states 
are not of single-particle nature and that the Green's 
function does not actually have any poles at uj = 0, but 
that there are zero-energy states that are collective exci- 
tations of many particles. In this case, the Green's func- 
tions will have zeros at the boundary to satisfy Eq. (11), 
but the actual system will still have zero-energy excita- 
tions. (However, we have not observed such a scenario in 
any of the examples considered in this paper). 



Regardless of the mechanism, the main conclusion is 
that if Ni is equal to 1, 2, or 3 modulo 4 (or if it changes 
by this amount across the boundary between two insula- 
tors), there must be zero-energy states of some kind at 
the edge of the system. 

III. ONE-DIMENSIONAL PEIERLS-HUBBARD 
MODEL 

Now let us consider a model of spin-1/2 fermions mov- 
ing in one dimension on a line with dimerized hoppings 
and on-site interactions, at half filling (one fermion per 
site). The Hamiltonian is given by 

'' ^(""' '') ("'• '-)• 
where (15) 

ft j> = 4<A CT ( 16 ) 

and Cj a creates a particle of spin er on site j — 1 . . . L. 
The ground state of this Hamiltonian is at half filling, 
that is, at one particle per site. This can be verified 
by observing that this Hamiltonian is invariant under a 
properly defined particle-hole transformation (see appro- 
priate detailed discussions in Ref. 14) 

c] a ^(-iyc^,c^^c]j-iy. (i7) 

Chiral transformations are particle-hole transformations 
combined with the time-reversal operation. The Hamil- 
tonian Eq. (15) is time-reversal invariant, so it is also 
chirally invariant. As a consequence of chiral symmetry, 
the Green's function satisfies Eq. (1). Chiral symmetry 
may still be broken spontaneously by the ground state, 
in which case the Green's function would then violate 
Eq. (1) despite the Hamiltonian's invariance. This ob- 
servation will be important later on. 

For what follows we need an explicit form of the gen- 
erator of the chiral transformation. One can check that 
it is given by 

£ = n (4 + (-^.t) (4 + i-^u) > (is) 

i 

which is both hermitian and unitary. It is straightforward 
to see that 

£ct CT £ = (-lyc^, ±c jt<7 £ = (-lyct^, (19) 
and that 

EH± = H, (20) 

which expresses the fact that the Hamiltonian Eq. (15) 
is chirally symmetic. 
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Chiral invariance together with time-reversal (and 
spin-rotation) invariance places this Hamiltonian in the 
symmetry class BDI, 1 although one should remark that 
normally the concept of symmetry classes is applied to 
nonintcracting systems only. But together with our ex- 
panded definition of topological invariants, we can extend 
this classification to interacting systems. 

In the absence of interactions, i.e., when U = 0, this 
Hamiltonian describes a particle hopping on the lat- 
tice with staggered hopping. Its topological invariant 
is straightforward to compute. Indeed, the Hamiltonian 
reduces to a matrix in this case, 



#=£^ 



jj' C j,cr C j' 



(21) 



To write this matrix in momentum space, we observe 
that the unit cell for this chain consists of two sites. By 
comparing Eq. (21) to Eq. (15) at U — we deduce that 



H(k) 



t + St + (t-St) e ik 

t + 5t+(t- St) e- ik 



(22) 

Here the lattice spacing is set to 1, and k £ (— tt, tt]. In 
the same basis, the matrix E takes its standard form 



E = 



1 

-1 



(23) 



By means of arguments presented earlier and given by 
Eqs. (7), (8) and (10), Ni simply measures the winding 
of the upper right corner 1 of % around the origin of the 
complex plane as k goes from — n to ir. This can be seen 
by substituting g = —'Hr x into Eq. (3) and observing 
that 



i dk d , 

N ^ 2l 2^dk h ^ 



(24) 



where z = t + St + (t — St) e lk . The factor of 2 in front of 
the integral signifies the fact that there are two identical 
copies of g, one for each spin. Finally, one immediately 
sees that 



m = 



2, for St < 
0, for St > 



(25) 



We have thus reproduced well-known facts about dimer- 
ized, noninteracting chains. In particular we see that 
they have a topological phase transition if St is tuned to 
be zero. If St is made space-dependent so that St > 
to the right of a given point in space and St < to the 
left of that point, there are zero-energy states localized 
in the vicinity of that point (as mentioned earlier, called 
Su-Schrieffcr-Hccgcr solitons 28 in this context). 

Now let us examine Hamiltonian (15) when U > 0. 
First take St = 0. In this case, this model is also very 
well understood. It is a Mott insulator, with a charge 
gap and gapless spin excitations described by a spin-1/2 
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FIG. 1. A fully dimerized finite Hubbard chain with \St\ = t 
consists of bonds of strength 2t (solid lines) and bonds of 
strength zero (not shown). The dashed box depicts the unit 
cell. In the case shown here, the topological invariant is and 
there are no boundary states. 



antiferromagnetic chain. 30 The strength of antiferromag- 
netic couplings between nearby spins is ~ t 2 /U if U ^S> t. 
However, St ^ dimerizes the spin-1/2 chain by mak- 
ing its bonds alternate in strength (while the charge sec- 
tor remains a Mott insulator). From the theory of spin- 
1/2 chains it is known that a perturbation consisting of 
nearest-neighbor spin-spin interactions with alternating 
sign is relevant. 31 From this discussion, it is natural to 
expect that even for U > the system is gapped for all St 
except St = 0, where a quantum phase transition occurs. 
This parallels the case without interactions. 

To understand the behavior of the topological invariant 
for U > 0, let us consider a limiting case when St = ±t. 
The chain then breaks up into disconnected clusters, each 
consisting of two sites only. The Green's function under 
these conditions can be found analytically by solving the 
two-site problem directly. However, it is not even neces- 
sary to solve the two-site problem to find the topological 
invariant. Indeed, the Green's function Gjji of this prob- 
lem is zero unless both j and f belong to the same two- 
site cluster. Its Fourier transform depends on whether a 
unit cell chosen previously when computing Eq. (22) co- 
incides with the cluster, or if in a given connected cluster 
one site belongs to one unit cell and the other site to an 
adjacent unit cell. The two cases are distinguished by the 
sign of St. In the former case, illustrated in Fig. 1, the 
Green's function is momentum- independent. Then g is 
also momentum-independent and N\ is obviously zero. In 
the latter case, illustrated in Fig. 2, the Green's function 
at zero frequency has the following structure in momen- 
tum space: 



512 e l 
-ik o 



(26) 



Indeed, we know that the matrix g must have zeros on 
its diagonal because g anticommutes with E, and it is 
easy to see that 512 must be a momentum-independent 









* 2t * * 







FIG. 2. A fully dimerized finite Hubbard chain as in Fig. 1, 
but with the alternate dimerization pattern; again, the unit 
cell is depicted by the dashed box. In the case shown here, 
the topological invariant is 2 and there are boundary states 
located at the isolated sites at both ends of the chain. 
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constant. Regardless of the actual value of 512, as long 
as the Green's function has neither poles nor zeros at 
uj = (that is, as long as gyi is neither infinity nor zero), 
the topological invariant in this case is N\ — 2, as can 
be established by substituting Eqs. (26) and (23) into 
Eq. (3). 

Therefore, we have established that Eq. (25) still ap- 
plies even when U > in the extreme case St = ±t. 
Earlier we saw that it also applies for all St if U = 0. 
Since we expect this model to be gapped, with the line 
St = in the U vs St phase diagram gapless, it is natural 
to conjecture that Eq. (25) applies for all St and all U. 
We will verify this in the next section. For now let us 
explore what knowledge of the value of the topological 
invariant implies for the boundary states of this insula- 
tor. 

According to the discussion in the previous section, we 
expect that a phase with N± — 2 has either zero energy 
states or zeros of the Green's function (or both) at its 
boundary. We would like to understand which of these 
possibilities is realized in our case. To do that, let us 
reexamine the case when St = ±i. 

When the unit cell consists of two sites at the ends 
of a nonzero bond, there are clearly no boundary states, 
as can be seen in Fig. 1 (and in that case, as we just 
saw, Ni = 0). In the other case, when Ni = 2, there 
will be a single unpaired site at each end of the chain 
(shown in Fig. 2). The Green's function for that site is 
straightforward to compute. The Hamiltonian of that 
single site is 

4o = U (4, t c, 0>t - (ct o4 c, o4 - , (27) 

where jo — 1 or j — L is the unpaired site at either 
end of the chain (L is the total length of the chain). The 
Hamiltonian acts in the space of four states — an empty 
site, a site occupied by a spin-up particle, a site with a 
spin-down particle, and a site filled with two particles — 
with energies U/A, —U/A, —U/A, and U/A, respectively. 
The ground state is then a site with either a spin-up or 
a spin-down fermion. These two states are degenerate. 
Therefore, it tells us that the ground state of the entire 
one-dimensional chain with two ends is four-fold degen- 
erate. These are what we can call the boundary states in 
the many-body context. 

Let us see how this is reflected in the Green's function. 
Choosing a state with either a spin-up or a spin-down 
particle as the ground state, we can calculate the Green's 
function directly from its spectral (Lehmann) decompo- 
sition. We find that the Green's function at this site is 

G io,CT;j ,T' = ( % "~v /2 1 ) • (28) 

Here the 2x2 structure of the Green's function is in 
spin space, with the first row and the first column of the 
matrix corresponding to the spin-up state if the ground 



state is chosen to be the state filled with the spin-down 
fermion, or vice versa. 

We see that the Green's function has neither poles nor 
zeros at ui = 0. The absence of poles is not surprising: the 
system is a Mott insulator and adding or removing a par- 
ticle should cost a finite energy, £7/2 in this case. Yet we 
need to see that this is compatible with, and indeed fol- 
lows, from the bulk-boundary correspondence, Eq. (11). 
At first glance, it contradicts the correspondence since we 
have a state with the topological invariant, N\ = 2, but 
without zeros or poles of the Green's function at ui = 0. 

Looking more closely, however, we observe that having 
a state with one particle on a site with a given spin vi- 
olates particle- hole symmetry. Indeed, the state ct Q ^ |0) 

goes into Cj g ^ |0) under the particle- hole transformation 
defined in Eq. (18), and vice versa; the doublet is sym- 
metric, not the states themselves. Mathematically, this 
is expressed in the fact that Eq. (28) does not satisfy 
Eq. (1). In other words, we observe that the presence of 
the many-body boundary states is reflected in the spon- 
taneous breaking of the particle-hole (and chiral) sym- 
metry. 

To summarize, we find that Ni = 2 does not result in 
either poles or zeros of the Green's function. This is not 
compatible with Eq. (11), as long as Eq. (1) holds. The 
way out is to recognize that the violation of Eq. (11) is 
only possible if there is a spontaneous breaking of the chi- 
ral symmetry and Eq. (1) is violated. In turn, this means 
that the ground state is multiply degenerate, which is a 
signature of the presence of zero-energy boundary states. 
This is how we can reconcile the value Ni = 2 with the 
absence of zero-energy, single-particle boundary states. 

But what if St is close but not quite equal to ±t and 
the chain is short enough that the boundary states can 
entangle and restore the chiral symmetry? Concentrating 
on the two ends of the chain with jo = 1 and jo = L with 
L the length of the chain, we find that the ground state 
must take the form 

(°4At - a * d lAi + ^lA.i + /^u^t) io) , 

(29) 

with arbitrary amplitudes a and (3 such that \a\ 2 + \/3\ 2 = 
1. Indeed, one can check that this state is an eigenstatc 
of E using Eq. (18), which for these two sites reduces to 

£= [ci ) t _a i.t][4 ) 4. _a u][cl,,t + Ci ) t][4,4. + Ci,4.]- (30) 

Here we ignore all the sites in the bulk, concentrating 
just on the two disconnected edge sites. 

Calculating the Green's function via its spectral de- 
composition yields 

G ]0 ,^> = <W ( ioj _ u/2 + lLJ + u/2 ) ■ (31) 

This function is indeed chirally invariant, satisfying 
Eq. (1). At the same time, it vanishes at uj = 0. 

We see that if the chiral symmetry is restored by con- 
sidering a finite chain with a unique ground state, then 
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the bulk-boundary relation, Eq. (11), is compatible with 
N% = 2 by way of the Green's functions having zeros at 
the edges of the chain. Even though this implies that an 
Ni = 2 chain does not have boundary states (something 
which would have been impossible in conventional topo- 
logical chains without interactions), this lack of boundary 
states is an effect of the finite length, and the boundary 
states are restored when the length of the chain is taken 
to infinity. 

Finally, while we have only looked at U > so far, a 
well-known mapping takes the half-filled Hubbard model 
at positive U to one at negative U . This is achieved by 
transforming 



(-l)^,t. c,- t (-1) 



(32) 



without changing C| and ct. Therefore, the whole discus- 
sion thus far applies equally well to the U < region of 
the phase diagram. 

Let us now summarize what we have learned by the 
direct analysis of the case St — ±t. If the topological in- 
variant is nonzero, one of three scenarios is realized. The 
first scenario is that the Green's function has poles at zero 
frequency, indicating zero-energy single-particle excita- 
tions. This occurs only at U = in our model. For U > 
the system is a Mott insulator, and adding a particle al- 
ways costs finite energy. The boundary states, if present, 
are particle-hole excitations, not single-particle ones; the 
single-particle Green's function cannot detect them and 
remains gapped. Nevertheless, either the Green's func- 
tion acquires zeros at the boundary or the chiral sym- 
metry is spontaneously broken by the boundary. These 
possibilities constitute the second and the third possible 
scenario, respectively. Either result signifies the presence 
of multiparticle, zero-energy boundary states in the limit 
of an infinite chain. 

Finally, we note that the ground state of the Peierls- 
Hubbard model is closely connected to that of a spin-one 
Heisenberg chain. To elucidate this point, let us add a 
ferromagnetic term to the Peierls-Hubbard Hamiltonian 
on every other bond. The combined Hamiltonian reads 



nu - 2 j + 



Here the spin S is defined by 

^ = O CjaCaa'Cjtj' , 



(33) 



(34) 



and cr is a vector of Pauli matrices. 

Since the Peierls-Hubbard Hamiltonian at large 
enough U can be thought of as a spin-1/2 chain, adding 



an explicit spin-spin interaction simply contributes to the 
dimerization of that spin chain. At the same time, it is 
clear that at very large negative J <C ~t 2 /U, spins will 
have a tendency to form a spin-1 moment out of two spins 
on sites 2j — 1 and 2j, converting the spin-1/2 antifcr- 
romagnetic chain into a spin-1 antifcrromagnetic chain. 
The spin-1 antiferromagnetic chain is known to be in the 
Haldane phase, which is gapped and has zero-energy edge 
excitations. 

If we switch off J while keeping St < 0, there is likely 
to be no phase transition, and the Haldane state will 
smoothly evolve into a state of a weakly dimerized spin- 
1/2 chain with topological invariant Ni = 2. This argu- 
ment supports the notion that the topological invariant 
in the Haldane state of a Heisenberg spin-1 chain is also 
Ni = 2. Thus, generically, we expect that the boundary 
states of spin-1 Heisenberg chains in the Haldane phase 
are due to the nonzero value of the topological invariant. 

We confirm this qualitative picture numerically in the 
next section. 



IV. DMRG ANALYSIS OF THE TOPOLOGICAL 
INVARIANT 

In this section, we discuss numerical results for 
the topological invariant in the Peierls-Hubbard model, 
Eq. (15), and for the extended model, Eq. (33), in 
which the ferromagnetic couplings lead to a Haldanc- 
like ground state. 32 We apply a Krylov variant 33 of the 
adaptive t-DMRG 19,20 to obtain the (real time) Green's 
function 



G(r,l) = (c l (r)cl /2 (O))0(r) 



(35) 



(here 9 is equal to 1 if its argument is positive and 
otherwise) for systems with up to L = 250 lattice sites. 
To keep the notation simple, we suppress the spin indices 
on the fermion creation and annihilation operators, but 
the Green's function in Eq. (35) is calculated for, say, 
spin-up fermions. Note that here and in the following 
we use the symbol r for the time in order to avoid con- 
fusion with the hopping strength t. To fix the units of 
energy and time, we set t = 1 throughout the section 
(although for clarity we keep t explicitly in some of the 
expressions below). We compute the topological invari- 
ant by obtaining the Fourier transform of this quantity, 
which we then use to analyze the expression in Eq. (3). 
However, since the relevant information is encapsulated 
in the phase of the Green's function, we avoid computing 
this quantity directly, which would require the computa- 
tion of the numerical derivative dkG(k,u), and instead 
analyze the winding of the chiral phase, which we define 
as 



V(k,u) = arg 



dre iUT E 

I even 



(36) 



Here, 1/2 effectively labels unit cells, k ranges from — tt to 
7r, and the formula is written assuming L/2 is odd. The 
winding V(k,uj) computed at uj = is nothing but the 
argument of z(k) introduced in Eq. (8). [More precisely, 
v{k) introduced there can be reduced to a number in our 
case, and V(k,0) is its argument.] 

Calculating V at lu — directly is not easy due to 
the slow convergence of the r integral, Eq. (36) . Instead, 
we use the analyticity and continuity of G to calculate 
its Fourier transform at a positive imaginary value of 
to, which ensures convergence. This does not introduce 
an error because the quantity we seek to calculate is an 
integer, and therefore is not sensitive to perturbations. 

Note that two complications limit the possible values 
of Imcj that can be used for this procedure. First, the 
magnitude of Imw must be such that, for the finite sys- 
tems under consideration, times after the perturbation 
has reached the edges are suppressed. Hence Imw > v/L, 
with v the speed of the fastest excitations in the system. 



(a) <n.>(T)forU = 2,5t = 0.3 
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FIG. 3. (Color online) (a) Time evolution of the local density 
(ni)(r) for a system with L = 22 sites, U/t = 2 and St = 
0.3. (b) Time evolution of the density on the central site 
for systems with L = 250 sites, St = 0.3, and for U/t = 
0, 0.3, 1, 2, 10. 



Second, Imw should not be too large compared to typi- 
cal energies of the system, since otherwise the structure 
of the Green's function will be determined by (the time 
cutoff provided by) Imu rather than by the eigenstates 
of the system. We conclude that Imu should be chosen 
to be of the order of the bandwidth, and, for the sake of 
simplicity, we choose u — i for all cases treated in this 
paper, without a noticeable effect on the value of the 
topological invariant. 

As explained below, finite-size effects are of minor im- 
portance. Thus, in order to reduce the computational 
effort, most of the results presented are obtained for sys- 
tems with only L — 22 lattice sites. For these computa- 
tions, we typically choose a time step dr = 0.05 and keep 
up to m = 500 density matrix eigenstates in the course of 
the time evolution, leading to a discarded weight of the 
order of 10 -5 at the end of the time evolution (r = 30). 
We adapt the number of basis states by choosing a dy- 
namical block state selection scheme (DBSS), 34 and fix 
the threshold for the quantum information loss (mea- 
sured by the Kholevo bound, see Ref. 35 for details) to 
X = 10~ 7 . Note, however, that, as we discussed above, 
the main properties of the chiral phase are determined by 
the behavior at short times (r < 5), so that the accuracy 
for these results is higher. 

Given a chiral phase, we simply plot V(k, i) as a func- 
tion of k and check if it winds around the unit circle as 
k goes over the Brillouin zone, from — 7r to tt. If it does 
not wind, then Ni = 0. If it winds once, Ni = 2. (The 
invariant is twice the winding because of spin.) 

In addition to computing the chiral phase, Eq. (36), 
we also analyze results for the real time evolution of the 
Green's function, Eq. (35), for the local particle density 
(n,i)(r), and for the local spin density (<Sf)(r), with a par- 
ticular focus on the role of edge states in the course of the 
time evolution. Due to the mapping (32), we treat only 
U > and expect the same behavior for negative U , with 
charge and spin interchanged. We start our discussion of 
the numerical results by considering the local observables 
(n,i)(r) and {S*)(t). In Fig. 3, we show a typical exam- 
ple for a system of L = 22 sites at U = 2, St = 0.3. As 
expected, at t = the density away from the center is 
very close to (rii) = 1 and does not show any signatures 
of end states for either positive or negative values of St. 
However, one realizes that due to the dimerization, the 
density distribution on the two central sites at the cen- 
ter of the system is unequal. In addition, as shown in 
Fig. 3(b), adding a particle leads to a larger density on 
that site with increasing U. For U = 10, near double 
occupancy is achieved. Interestingly, in the time evolu- 
tion of the density at the central site [Fig. 3(b)] the curves 
seem to intersect at r ~ 0.5 for values of U small enough. 

Edge states appear in the time evolution of the local 
spin density (Sf) (r). As can be seen in Fig. 4, for positive 
values of St, we find that (S?)(t = 0)k0 away from the 
central sites. However, for negative values of St, bound- 
ary states appear. For small values of U, Fig. 4(b), no 
signature of boundary states can be seen at the beginning 
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of the time evolution. However, in the course of the time 
evolution, local spins (pointing in opposite directions) 
seem to emerge at the two boundaries. For larger values 
of U [Fig. 4(c)], these end spins are already visible in 
the initial state and appear to be close to fully polarized, 



(a) < Sf >(t) for U = 2, 5t = 0.3 




(b) < Sf >(t) for U = 0, St = -0.3 
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site i 



(c) <Sf >(x) for U = 2, 5t = -0.3 
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FIG. 4. Time evolution of the local spin density (5 , f)(r) for 
(a) U = 2 and St = 0.3, (b) U = 0, St = -0.3, and (c) 
U = 2, St = -0.3. 



again with opposite orientations at the two boundaries. 
This is consistent with the picture of boundary states 
for St < 0, a picture we subsequently confirm with the 
computation of the topological invariant. 

Now we turn to the time evolution of the Green's func- 
tion. In Fig. 5, we show typical results for the case of a 
system with L = 22 sites, U = 2, and St — 0.3 for times 
up to t = 30, the maximum time reached. Again, the 
dimerization leads to an asymmetry in the initial state. 
The perturbation in the center spreads through the sys- 
tem with a typical velocity that depends on the values of 
U and St. No signature of boundary states can be seen 
for all values of U and St treated. As can be seen, at 
a time t 10, the perturbation reaches the boundary 
and gets reflected. This can lead to significant finite-size 
effects in the Fourier transform at uj = 0, giving another 
reason why one should set u> = i, as discussed earlier. 
This effectively suppresses contributions to V from late 
times, which are influenced by reflection from the bound- 
aries, and so better captures the behavior of the infinite 
system. We find that for L = 22, all values of U < 10, 
and —t<St< t, the perturbation reaches the boundary 
at a time r w 5 or later. By choosing uj = i, which is a 
value of the same order of magnitude as the bandwidth, 
results for G(r, I) at times r > 5 have a very small weight 
(< 1%) in the Fourier transform Eq. (36), so that finite- 
size effects should become minimal, as further discussed 
below. 

Based on this, we present in Fig. 6 results for V(k, uj — 
i) for St = ±0.3 and U = 0, 2, 10. For positive values 
of St, we do not expect boundary states to be present, 
and there should be no winding of V(k, i), while for neg- 
ative values of St, one winding of V{k, i) should be ob- 
tained. This is indeed the case: as can be seen in Fig. 6, 
for positive as well as negative values of St the function 
V(k,i) appears to be periodic. As seen in Fig. 6(a), 



G(x,l) for U = 2, St = 0.3 




FIG. 5. (Color online) Time evolution of the Green's function 
for a system with L — 22 sites with U = 2, St = 0.3. In the 
case shown, up to m = 750 density matrix eigenstates were 
kept during the time evolution. 
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FIG. 6. (Color online) Chiral winding V{k, i) [Eq. (36)] for 
systems with L — 22 sites (large symbols, solid lines) and 
L = 250 sites (small symbols), (a) Results for St — 0.3 and 
U = 0, 2, 10. (b) Results for St — —0.3 and the same values 
off/. 



for St > 0, the values covered are restricted to a region 
7r/2 < V(k, i) < 7r, so that the winding number is zero. 
For negative values of St, however, the values cover the 
full range from — it to n, so that the winding is equal to 
1. Thus, we arrive at the main result of this paper: even 
in the presence of interactions, the topological invariant 
is equal to zero or 2 (one winding per spin) , and its value 
reflects the presence of boundary states. 

Now we turn to finite-size and finite-time effects. In 
Fig. 6, we compare our results for V(k, i) for U = 0, St = 
±0.3 to results with L = 250 sites. For L — 22, times 
t = 15 or larger were reached, so that, according to the 
above discussion, finite-time effects should not be a major 
issue. For L = 250, the computations are more demand- 
ing and only times of the order of r ss 3 were reached 
in the cases shown. Apparently, on the scale of the plot, 
finite-size and finite-time effects seem to be absent. Wc 
further analyze this in Fig. 7 in which we compare results 
for system sizes L = 22, 50, and L = 250 for U = 10 and 
St = 0.3. For L = 22 and L = 50, times of at least 
t = 12.5 were reached, while for L = 250, only times 
of t = 1.6 were reached. By comparing the results for 
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FIG. 7. (Color online) Chiral winding V(k,i) at U = 10 and 
St — 0.3 for systems with L — 22, L = 50 and L — 250 lattice 
sites. While times r = 12.5 or larger were reached for the 
smaller systems, only times r = 1.6 were reached for L — 250 
sites. 



L = 22 and L = 50, we see that finite size effects seem 
to be practically absent also in this strongly interacting 
case. However, the results for L = 250 seem to be shifted. 
We associate this with the small times that were reached. 
Importantly, this finite-time effect does not affect the 
overall behavior, and the winding number is correctly 
obtained. We therefore conclude that computing the chi- 
ral phase is, at least in the present case, a very stable 
numerical procedure and can be performed with a rather 
moderate numerical effort for rather small systems and 
short times. Finally, we discuss results for fixed values of 
U and varying St, shown in Fig. 8 for U = 10 and L = 22. 
As can be seen, for all values of St < —0.05 shown, one 
winding is present, while for all values of St > there 
is no winding. This is in agreement with the well-known 
phase diagram of the Peierls-Hubbard model, 26 in which 
the spin gap closes on the line (U,St = 0), signifying a 
phase transition. Via computing the topological invari- 
ant, Fig. 8 shows that this phase transition connects a 
topologically trivial phase for St > to one with bound- 
ary states at St < 0. The phase diagram of the Peierls- 
Hubbard model with the two phases identified by the 
value of the topological invariant is shown in Fig. 9. 

Based on these insights obtained in the Peierls- 
Hubbard model, we can now attempt to apply a similar 
procedure to other systems. Here, we discuss the chiral 
phase for the extended model (33) with spin-exchange in- 
teractions. At U > 0, the chain is a Mott insulator and, 
just as in the discussion of the model (15), all its physics 
is in the interaction between the on-site spins. For suf- 
ficiently large ferromagnetic coupling J/t < —t 2 /U, wc 
expect the formation of effective spin-1 objects, leading 
to a ground state similar to that of a spin-1 Heisenberg 
chain with boundary states. In contrast, antiferromag- 
netic coupling enforces the formation of singlets, so that 
the situation should be similar to the above discussion 
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for the Peierls-Hubbard model at St > 0. In Fig. 10, we 
present results for V(k, i) in the strongly interacting case 
at U = 10 with L = 22 lattice sites for a ferromagnetic 
J = — 2 and an antiferromagnetic J = 2 with St = 0. We 
observe a winding of 1 (corresponding, as before, to the 
topological invariant Ni — 2) if J = —2 and Ni = if 
J = 2. (Note that i 2 /[7 is 1/10 and is much smaller than 
J.) This confirms that for negative J, where we expect 
this system to be a spin-1 Heisenberg chain, boundary 
states are present. 
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We could now investigate other aspects of this sys- 
tem, such as the critical value of J at which the invariant 
changes and the system goes through a quantum phase 
transition or the dependence of that transition on St. 
However, since this would lead us too far afield, we leave 
the investigation of the full phase diagram of the model 
Eq. (33) to future studies. 



FIG. 9. Phase diagram of the Peierls-Hubbard model given 
by Eq. (15), with the value of the topological invariant Ni 
indicated. The line St = separates two gapped phases, 
topological and non-topological. The line (7 = describes 
the usual noninteracting ID topological insulator. 



ANALYSIS OF THE FIDKOWSKI-KITAEV 
MODEL 
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FIG. 8. (Color online) Chiral winding V(k, i) for systems with 
L = 22 sites, U — 10 and different values of St: (a) St > 0; 
(b) St < 0. 



As discussed in the introduction, a one-dimensional 
interacting system with the topological invariant that is 
a multiple of 4 may have no boundary states whatso- 
ever. This means that two Hamiltonians that are topo- 
logically distinct at the quadratic, noninteracting level 
can be adiabatically connected by adding an appropri- 
ate interaction. Here we demonstrate that such a system 
must have Green's-function zeros at its boundary. (Cor- 
respondingly, zeros develop at zero energy somewhere 
along the path in the parameter space that deforms one 
such system into another one with a distinct topological 
invariant.) 

To carry out the demonstration, we present an analysis 
of the Fidkowski-Kitaev model, 23 which consists of two 
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FIG. 10. (Color online) Chiral winding V(k, i) for the ex- 
tended model with spin-exchange terms, Eq. (33), for systems 
with L = 22 sites, U = 10, and J = 2 (antiferromagnetic case, 
red squares) and J = —2 (ferromagnetic case, blue circles). 
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FIG. 11. Graphical representation of the Fidkowski-Kitaev 
model of Ref. 23: Two dimerized Hubbard chains are coupled 
via a "rung" Heisenberg exchange term. 



Peierls-Hubbard chains coupled by a Heisenberg interac- 
tion. The topological invariant of this model is either 
or 4, even though it has only one phase. The Fidkowski- 
Kitaev Hamiltonian has a large symmetry that is more 
apparent when it is rewritten in terms of real (Majorana) 
fermion modes, which will be largely hidden in the treat- 
ment that follows. 

The Heisenberg interaction between the two chains a 
and 6, described by fermion operators and is 



JS a ■ S;,, 



(37) 



where the spin operators are defined in terms of the 
fermion operators as before. We assume J > (i.e., 
is antiferromagnetic) in the following, so that "rung sin- 
glet" is energetically favorable. The full Hamiltonian is 
depicted graphically in Figure 11. 

The addition of the Heisenberg term allows a smooth 
interpolation between the two topologically distinct 
phases with signJt = ±signi. Essentially, it provides 
an extra source of dimerization in the system that is lost 
when St = 0, thereby keeping the spin waves gapped 
(in the case U > 0). Let us outline the argument of 
Fidkowski and Kitaev to this effect and then reconcile 
this with the topological invariants of the single-particle 
Green's function. 

As before, the simplest limiting case of the system is 
St = ±t. In this case, the topologically nontrivial phase 
for J = has zero-energy states located purely at the 
ends of the chain, with no tunneling. As we have seen, 
this is true even for nonzero U, when a single-particle 
gap opens but two degenerate ground states per chain, 
with a spin-1/2 at the end and no dynamics, remain. The 
Heisenberg term then causes a singlet to form from these 
spins, with a gap to the triplet excitations. 

The bulk part of the fully dimerized chain is gapped, 
and adding the Heisenberg interchain interaction does 
not change much qualitatively. After all, the dimers 
of the Hubbard chain can effectively be described by a 
Heisenberg interaction with J' ~ t 2 /U. Adding the ex- 
tra Heisenberg term allows the dimerization to interpo- 
late smoothly from the chains to the rungs. This means 
that, for J/t > and U/t > 0, the tunneling can be 
turned off entirely without closing the gap, making the 
adiabatic continuation between the two noninteracting 
phases, which is confirmed by the lack of zero-energy 
modes at the ends, possible. 

The question is how can this be consistent with the 
fact that the noninteracting phases and the interacting 



phases with J = are distinguished by a topological 
invariant? As before, the answer is that the bulk single- 
particle Green's function must develop zero- frequency ze- 
ros when St = 0. Similarly, the Green's function for the 
end sites must have zero-frequency zeros for St < 0. This 
situation is distinct from that of a single Hubbard chain, 
in which the Green's function breaks chiral symmetry due 
to the degenerate ground state; here the ground state is 
unique. 

To compute the Green's function, we must determine 
the ground state, the spectrum, and the matrix elements. 
We will only consider the Hamiltonian i?Hub + H spnl ; it 
describes the end state of the chain with St = —t or the 
bulk for t = St = 0. When J > the ground state, with 
energy E gs = —3 J/4 — U, is 
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(38) 



where |0) is the Fock vacuum of the a and b fermions. 
The single-particle excitations above this ground state, 
given by acting on the ground state with a single creation 
or annihilation operator, are all degenerate, with energy 
E\ = —U/2; there are eight such states. The remaining 
states have energies J/4 — U (triplet) and 0. Tuning to 
the point J/4 — U increases the symmetry of the model, 
as discussed at length by Fidkowski and Kitaev. 

The single-particle Green's function of the decoupled 
chain is given by 



Gij{u) = (gs\ fi(iuj + E g 
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(39) 



Here the fermion operator /, takes on the four values a-j-, 
aj,, bf, and b±. Given the properties listed above, the 
Green's function evaluates immediately to 



Gij(u) 
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ius - 3 J/4 - U/2 iuj + 3 J/4 + U/2 



(40) 

This function satisfies chiral symmetry [Eq. (1)], but does 
not have a pole at zero frequency because the single- 
particle excitations are gapped. Instead, it has a zero, as 
Gij(Q) = 0. This confirms the reasoning outlined above. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper, we have demonstrated the utility of cal- 
culating the topological invariant for interacting topolog- 
ical gapped systems, working with the example of spinful 
fermions hopping on a one-dimensional lattice. While the 
invariant is no longer directly related to conductivity or 
other responses of the system, it can still be used to de- 
duce whether zero-energy boundary states are present, 
thanks to the bulk-boundary correspondence, Eq. (11). 
The invariant can be computed numerically; in this pa- 
per we accomplish this with the DMRG method. One 
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advantage of computing topological quantities numeri- 
cally is that, as integers, they are not strongly prone to 
numerical errors. 

One could, in principle, ask whether direct numeri- 
cal evaluation of the boundary states is no more difficult 
than evaluating the topological invariant, or whether it 
may even be advantageous. We would like to point out 
that direct evaluation of this sort is more prone to nu- 
merical errors. The boundary states are susceptible to 
finite-size effects. While they could be lifted away from 
zero numerically, determining whether they are topolog- 
ically protected might then not be easy. The topologi- 
cal invariant is robust and is only weakly susceptible to 
numerical errors. If it is found to be nonzero, the zero- 
energy boundary states are guaranteed to exist in the 
large-size limit. 

While direct numerical determination of the topologi- 
cal invariant such as Eq. (3) can be problematic due to 
the numerical errors associated with integrating deriva- 
tives of functions determined numerically, in the present 
one-dimensional case this was not necessary. Instead, 
we evaluated the winding associated with the topological 
invariants by inspection. For completeness, let us note 
that, if desired, we could have evaluated the invariant 
without having to inspect the graph of the chiral phase, 
Eq. (36), visually. Instead, we could have found all the 
solutions ki of the equation 



of a matrix G{k, ui) given by 



V{ki,0) = V , 



(41) 



where Vo is an arbitrarily chosen number between — tt and 
7r. Given the set of ki which solves this, we can compute 



dV(k,0) 



Ok 



k=kt 



(42) 



24tt 2 ? 



t aM tr J dud 2 kG- 1 d a GG- 1 d p GG- 1 d 1 G. 



(43) 

Here a, j3, and 7 are summed over uj, k x , and k y . 
This is equivalent to computing the Chern number 
if there are no interactions 37 and can be reduced to 
a two-dimensional Berry-curvature integral even with 
interactions. 12,13 Evaluating the derivatives and the in- 
tegrals in Eq. (43) numerically is problematic. Instead, 
in the important case where G is a 2 x 2 matrix, one can 
parametrize it by writing it as a sum over a unit ma- 
trix and the three Pauli matrices with coefficients vq, v%, 
1)2, V3. Since the overall normalization is irrelevant, this 
corresponds to three parameters (61,02,63) that are 
functions of (ui,k x ,k y ). Then, given G, we could evalu- 
ate the Jacobian at some special value of 0, an analog 
of Vo in Eq. (41). The sum of the signs of the Jacobians 
so computed is equal to Eq. (43). This method works if 
G is a 2 x 2 matrix. Its generalization to the case where 
G is a larger matrix is not known to us, but should, in 
principle, exist. 

Now that we have established that the method of 
topological invariants is useful for studying interacting 
fermionic systems in one-dimensional space, it would be 
interesting to apply it to other interacting topological 
insulators. One possible direction of further research 
would be to study two- and three-dimensional topolog- 
ical interacting systems. It might also be interesting 
to further apply this method to other one-dimensional 
problems, for example, spin chains and ladders accessi- 
ble to the DMRG. It also would be worthwhile to clar- 
ify the relationship between this method and recently 
discussed symmetry-protected topological orders in one- 
dimensional space. 38 



This works for almost all Vo and is Vb-independcnt. (It 
fails for those Vq for which V(k, 0) has a vanishing deriva- 
tive at k being equal to one of the ki.) If the derivative 
is too small to determine its sign dependably, a different 
Vo can be chosen. As elsewhere throughout the paper, 
the prefactor 2 has to do with the two spin components 
of spinful fermions. 

Importantly, the formula Eq. (42) has a natural coun- 
terpart in higher dimensions. 36 For example, in two spa- 
tial dimensions one might want to evaluate the winding 
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